Density profile of a harmonically trapped ideal Fermi gas in arbitrary dimension 



o 
o 

(N 
>v 



Erich J. Mueller 

Laboratory of Atomic and Solid State Physics, Cornell University, Ithaca, New York 14853 

(Dated: May 21, 2004) 

Closed form analytic expressions are derived for the density profile of a harmonically trapped 
noninteracting Fermi gas in d dimensions. Shell structure effects are included to leading order in 
1/N, where N is the number of particles. These corrections to the local density approximation scale 
as 8n/n ~ N~ a , where a = (1 + l/d)/2. 

PACS numbers: 03.75.Ss 
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Experiments on quantum degenerate Fermi atoms have 
motivated a series of theoretical studies of the basic prop- 
erties of zero 0, Q, S S S an d finite Q tempera- 
ture noninteracting fermions in inhomogeneous poten- 
tials. Here we derive closed form analytic expressions 
for the ground state density of a harmonically trapped 
Fermi gas in d dimensions. 

Theoretical studies of noninteracting particles are cru- 
cial for understanding experiments on nonresonant Fermi 
gases @- A typical length scale for interactions is 
ro ~nm, and a typical interparticle spacing is n" 1 / 3 ss 
/im3> ?"o- Thus, even in multiple component gases, where 
s-wave collisions are allowed, the interaction energy per 
particle E- la t « h 2 rQn/m is small compared with the 
Fermi energy Ef w Ti 2 n 2 / 3 jra. Therefore the interac- 
tions can be ignored for calculating gross features of the 
ground state structure. This separation of energy scales 
is even more dramatic in a spin-polarized gas, where s- 
wave collisions are forbidden, and therefore p-wave colli- 
sions dominate. The cross-section for p-wave collisions is 
down by a factor of (fc/ro) 4 relative to s-wave. 

Although our calculations are for arbitrary dimension, 
d, a particularly important case is d = 1 where theoretical 
predictions about noninteracting fermions can be applied 
to hard-core bosons @. 

Most experiments are performed in harmonic traps, so 
it is natural to consider the density profile of a gas of 
N fermions in a potential V(r) = muj 2 r 2 /2. Since we 
ignore interactions, each spin component is independent, 
and it suffices to consider the spinless (or spin-polarized) 
case. The simplest expression for this density profile 
comes from the local density (Thomas-Fermi) approxi- 
mation, where the trap potential gives a spatially depen- 
dent chemical potential /i(r) = k 2 {r)/2m = ^ — V{ r ), 
and the local density is derived from the relationship be- 
tween density and chemical potential in a (/-dimensional 
homogeneous system, n = (kf /2*Jir) d /T(d/2 + l), result- 
ing in a density profile, 
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where R 2 JZ 2 = 2{d\Nf/ d is the radius of the cloud,and 
I 2 = h 2 /mtu 2 is the oscillator length. Although this is 



an excellent approximation to the true ground-state den- 
sity profile, it misses "shell effects" resulting from the 
statistical correlations in an ideal Fermi gas. These shell 
effects give rise to corregations on top of this smooth pro- 
file which are analogous to the Friedel oscillations which 
occur in the density of a uniform Fermi gas near an im- 
purity [ljj. As we show below, these corregations scale 
as Sn/n ~ N~ a where a = (1 + l/d)/2. These devi- 
ations are therefore only significant for small numbers 
of par ticles. Despite their small size, previous estimates 
yj, [llj suggest that these deviations are experimentally 
observable. 

Several authors have quantified these density fluctua- 
tions. In one dimension, Husimi |l2| derived a simple 
closed-form expression for the exact density profile in 
terms of Hermite polynomials. Further discussion and 
applications of this approach can be found in Us- 
ing a completely independent formalism, Vignolo et al. 
[l| recently derived an efficient numerical technique for 
calculating this profile. Vignolo and Minguzzi later gen- 
eralized their numerical approach to higher dimension . 
Also in higher dimensions, Brack and Zyl Q derived an 
analytic expression for the exact density profile in terms 
of a sum over Laguerre polynomials. Brute force numer- 
ical calculations of three dimensional shell effects can be 
found in |l4| . including some discussion of corrections 
due to interactions. 

Unlike these prior works, we will not attempt to de- 
rive an exact expression for the density profile, rather 
we will produce an asymptotic series in powers of 1/N. 
Consequently we are able to derive simple analytic for- 
mulae. Near the center of the trap, our results agree 
with those derived by Gleisberg et al and Brack and 
Murthy [|J. Our expressions match the exact density 
profile throughout the cloud, even for surprisingly small 
numbers of particles: N ~ 1. 

Density Generating Function: The single parti- 
cle states of a d dimensional harmonic oscillator are 
described by the quantum numbers {ni, . . . , rid}, cor- 
responding to a real-space wavefunction of the form 
ip(xi,...,Xd) = (f) ni (xi)---(f> nd (x n ), where <pj(x) = 
H J (y)e"y 2 / 2 /(2^jUy/^) 1 / 2 is the eigenfunction of the 
one-dimensional harmonic oscillator, the dimensionlcss 
length is y = x/t, and Hj(y) is the j'th Hermite 
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The generating function for the contribution to the den- 
sity in d dimensions from the M'th shell is the product 

d 

G(r,t) ="£* d m (r)t™ = Y[W( Xj ,t). (6) 

m j=l 



FIG. 1: (color online) Contours of steepest descents in the 
complex t plane for regions A,B,C,D (see text). Poles (at t = 
0, 1) are represented by red squares, saddle points by yellow 
discs, essential singularities (at t = — 1) by blue triangles, and 
branch cuts (extending to ±oo from t = ±1) by green wiggly 
lines. The dashed circle corresponds to \t\ = 1. 



polynomial. This single-particle state has energy E — 
Tiu>Y2j( n j + 1/2)- The ground state of N fermions in 
a harmonic trap is only unique if an integral number of 
shells are filled, requiring 
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where M = 0,1,2,... indexes the shell. The Thomas- 
Fermi radius is then R 2 ti /t 2 = 2M+(l+d)+0(l/M). The 



exact density will be denoted p^f ( r ) = El'=i ^'W 



where o£?(r) = E„ 1+ ..+„ d = m 1 0m i x i) ' ' -<t>n d {xd)\ 2 is 
the density contribution from the m'th shell. In the re- 
mainder of this paper, we scale all lengths by t. 

We will extract the large M structure of the density 
from the generating function, G(r, t) — YIm Pm ( r )t M ■ 
We construct G by noting that the Hermite polynomials 
obey 



F(x,z) 
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exp(-z 2 + 2zx). (3) 



To construct expressions involving the density contribu- 
tion from each orbital, |(/>j(a;)| 2 , we need a generating 
function for Hj{x) 2 . We therefore introduce 



W(x,s) = \du\dve- {u+v )/s F(x,u + iv)F{x,u-iv) 
H 3 {x) 2 s 2 i+ 2 



(4) 
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vT - 4s 2 
A simple rescaling then gives 

w{ X ,t) = Y.\^ x )\ 2t 



3 n 3 /H 



e~ x W{x,s/t/2). (5) 



The density for M filled shells then obey 

G(r,t) 

M 
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exp 
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which shows (as would be expected from symmetry) that 
the density is only a function of the distance r = |r|. 
Asymptotic Expansion: The relationship J7J is in- 
verted by performing a contour integral around the origin 
in the complex t plane, 



P M \r) 



dt 
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Deforming the contour to follow the path of steepest 
descents, the large M asymptotics are calculated from 
the portion of the integral which passes over saddle 
points. The path of steepest descents is illustrated in 
figure n f° r four different regimes: A : r 2 /i? 2 f < -R^ 4 , 
B : i?" f 4 < r 2 /R 2 ti < 1 - R^'\ C : \r 2 / R 2 ti -l\< R^'\ 
and D : r 2 /R 2 t — 1 > i? tf 4 ^ 3 . These regimes are distin- 



guished by the number and locations of the saddle points 
which are found on the steepest descents contour. 

In all cases the integral is dominated by a single sad- 
dle point on the real axis at t = to with < t n < 1. For 
r/Rtt < 1 — C(i? 4 / 3 ), this saddle point approaches the 
pole at t = 1 as to = 1 + 0(1/M), and one cannot use 
the standard saddle-point formulae to evaluate its con- 
tribution. However, straightforward analysis shows that 
for a sufficiently well behaved 0] function f(t), one can 
approximate 
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plus terms of relative order 1/M. The integral is taken 
on a contour from t = ~ioo to t = +ioo with Re(t) < 1. 
Near t — 1, G{r,t)/t M ~ x has the same form as the inte- 
grand in ©, with M/(l) = (d/2) log(27r) and M/'(l) 



2 )/2 + (M° ) , where R 2 tl = 2M + 1 + d. Thus the 



leading order behavior coincides with the Thomas-Fermi 
result in equation fljl. 

The other saddle points are found by setting the log- 
arithmic derivative of the integrand in (JSJ equal to zero. 
The resulting cubic equations yields two additional sad- 
dle points, denoted t±. In region A, t± — ► — 1; in B, t± 
are isolated saddle points with \t±\ « 1; in C all three 
saddle points converge at t = 1, and the contour of steep- 
est descents no longer passes through t±; and in D, all 
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FIG. 2: Density profile of one-dimensional noninteracting 
fermions in a harmonic trap with frequency w. Profiles with 
N=l to N=ll particles are shown. Density p and position r 
are measured in units of the trap length I — yjh /mui. Thick 
gray lines: exact result; thin black lines: Eq. (I1H . The exact 
and approximate results lie on top of one another except near 
the edge of the cloud. 



three saddle points lie on the real axis, and to moves away 
from t = 1, becoming an isolated saddle point. In the last 
two regions we must explicitly consider the corrections to 
© ■ Note that most of the cloud lies within region B and 
we therefore first consider the Friedel oscillations in this 
regime. We then move on to regions A, D, and C. This 
order is chosen so that we can sequentially use previous 
results. 

Region B. Integrals passing over isolated saddle points 
can be evaluated through the standard formula, 



dse M ^ = 



2tt 



-Mf"(s) 



+ 0(1/M 1 / 2 )), 



(10) 

where s is the saddle point, characterized by f'(s) = 0, 
and the integral on the left represents the contribution to 
the integral along the path of steepest descents near the 
s. To simplify the algebra we change variables, letting 



t 



Setting the logarithmic derivative of G equal 
to zero, one finds that the saddle points t± correspond 
to s± = ± arccos(r/i? t f ) + 0(l/M). In this new vari- 
able, the integral in (|SJ) maps onto (|10|) with M f(s) = 
ir 2 tan s - isR 2 { - (d/2 + 1) log(7t) - (d/2) log(2 cos(s)) - 
(l + d/2) log(2i sin(s)). Adding the contributions of these 
two saddle points, one arrives at out central result: away 
from the center and edges of the trap, the leading order 
form of the density oscillations is 



A, 



-(3 + d) 

(l - f 2 ) 4 cos(i?. 2 
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0(l/M)). 



M (ii) 

(12) 
(13) 
(14) 
(15) 



By examining Eq. (|11|1 . one sees that Sp scales as l/\/~M 
for large M, while the Thomas-Fermi density scales as 
M d / 2 . The accuracy of (|llfl is illustrated by the examples 
shown in figure |21 



Region A. When r — > the poles t± approach one an- 
other at t — — 1, and no longer contribute as isolated 
saddle points. The width of each saddle point scales as 
(As) 2 ~ Rtf/ r as r ~~ * 0, while the distance between the 
saddle points (s+ — s_ ~ arccos(r/i? t f)) approaches zero 
as r/Rtf. Therefore the isolated saddle point approxima- 
tion breaks down when r 2 < R^ { 2 - 

To analyze this limit, we expand G for small r and 
U = t + 1. The leading order behavior for 5p is 



8p / 
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with w = M + 3/2 + d/4 w i? 2 f /2. By deforming the 
contour, this integral may written in terms of a Bessel 
function, J t _ l (z) = (2iri)~ 1 J duu - ^ -1 exp(z(u — l/u)/2), 
to yield, 
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(-1) 



A/ 



2 d ir d / 2 



-l-d/2 



J d/2 ^(2fR 2 { ). 



(17) 



This result was previously found in d dimensions by 
Brack and Murthy |3J, and in 1 dimension by Gleisberg 
et al. : .5j- By considering the long distance asymptotic 
behavior of the Bessel function J v {z) ~ \j2j-KZ co&(z — 
■nv/2 - 7r/4)(l + 0{z^ 2 )), one sees that for R~ 2 «f« 

2/3 

i? tf the solutions bp a and Sps coincide. 
Region D. As one makes r larger, the saddle points 
which contribute to Sp follow the circle \t\ — 1, and 
as r ^ i? t f they approach the pole at t — 1. For 

— 1/3 

r > Ra + 0(R t{ ) all three saddle points move onto 
the real axis; two with t > 1 and one with Q < t < 1. 
The path of steepest descents only passes over the last 
of the saddle points. The Thomas- Fermi approximation 
has no applicability here, so we calculate the full density, 
rather than the deviation from Thomas-Fermi. Since the 
relevant saddle point is isolated, we can use the standard 
formula from equation l|10|l . Similar to our approach in 
region B, we change variables, setting t = e 2w , so that 
the density is 



Pd{t) 



1 



i2 1+d ir 1+d / 2 



dw e , 



(18) 



where 9 = —R 2 { w + r 2 tanhu> — (1 + d/2) logsinh(— w) — 
(d/2) log cosh w. To leading order, the saddle point is at 
Wo = — arccosh(f), yielding a gaussian decay, 



(1 



Pd = A D r~ 
A D = n -(d+i)/2 2 -(d+i/2) R -^ 

arccosh(f). 



A 
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exp(-i? 2 f A) (19) 
(20) 
(21) 



Comparing the width of the saddle point with the dis- 
tance between wo and the pole at w — 0, one finds that 
this expression is valid when f ^ 1 + R 



-4/3 
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The density profile in this region can also be estimated 
from a modified Gross-Pittaevskii type functional |l6| . 
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FIG. 3: Density profile of four filled shells (M=3, N=10) 
of non-interacting fermions in a two dimensional harmonic 
trap. Top and bottom panels respectively show density p and 
deviation 8p from a Thomas-Fermi profile. For clarity, the 
various approximate solutions are shown only on the left or 
right half of the graph. Thick gray line: exact; solid black 
line (right): Eq. JTTJ; dotted line (right): Eq. dashed 
line (left): Eq. (H); dot-dashed line (left): Eq. fig) . 

However, such an a ppr oach can lead to unphysical inter- 
ference predictions |17| . 

Region C. Between regions B andZ? (1-1/ R 4 J 3 < f 2 < 

1 + all three saddle points are in close proximity 

to t = 1 (w — in the variables used to discuss region D). 
Linearizing 9 about w — 0, one finds to leading order, 

6 = -(r 2 - R 2 { )w - (1 + d/2) log(w) - r 2 w 3 /3. (22) 

By introducing a small parameter, r\ = (3 1 ^ 3 r~ 2 / 3 (i? 2 f — 
r 2 )), and changing integration variables to v = r 2 w 1 /3, 
we can write the integrand in a power series in r\. Each 
term can be integrated by using a variant of equation 
yielding 

« - ^ ,a i: ro+1 , r( ^ +1 _ J/3) p3> 

A c = 3 -(i+<i/6) 2 -rf 7r -rf/2 (24) 

This expansion is rapidly convergent when r\ < 1, corre- 
sponding to 1 - l/i?t/ 3 < r 2 < 1 + l/i? t 4 / 3 - 
Summary: We have analytically calculated the density 
profile for N noninteracting spin-polarized fermions in a 
e?-dimensional harmonic trap. Four different expressions 
are needed to describe all parts of the cloud. The results 
from all four regimes are all sketched in figure [3] for the 



case of ten particles in two dimensions (corresponding to 
four filled shells). 
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